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Abstract. - The lattice Boltzmann algorithm efficiently simulates the Navier Stokes equation 
of isothermal fluid flow, but ignores thermal fluctuations of the fluid, important in mesoscopic 
flows. We show how to adapt the algorithm to include noise, satisfying a fluctuation-dissipation 
theorem (FDT) directly at lattice level: this gives correct fluctuations for mass and momentum 
densities, and for stresses, at all wavevectors k. Unlike previous work, which recovers FDT only 
as fc 0, our algorithm offers full statistical mechanical consistency in mesoscale simulations 
of, e.g., fluctuating colloidal hydrodynamics. 



The lattice Boltzmann equation (LBE) is a widely used lattice formulation of fluid me- 
chanics [1]. It offers a faithful discretization of the Navier Stokes equation of isothermal, 
incompressible fluid flow, and is very well adapted to parallel computation [2]. While used 
for large-scale fluid dynamics simulations such as flows around aircraft [3] , the LBE approach 
is particularly adapted to simulating mesoscopic problems [4]. These include, e.g., porous 
medium flows and flows of complex and multicomponent fluids with microstructure [5-8]. 
The latter can be modelled using various extensions of the basic algorithm for a single com- 
ponent fluid as considered here [8-10]. 

However, the Navier Stokes equation, and with it the LBE, ignores thermal fluctuations. 
While these may safely be ignored in macroscopic fluid-dynamical flows, at mesoscopic length 
scales they form an essential part of the physics [11]. This applies even in linear problems 
such as the Brownian motion of a colloidal particle suspended in a simple fluid: if that fluid 
is simulated using the LBE, no Brownian motion occurs [8]. Fluctuations are also central 
to nonlinear phenomena such as mode-coupling effects and long-time tails [12]. By the same 
token, extensions of the LBE to fluid mixtures [9] and amphiphilic solutions [10] cannot address 
critical phenomena, where fluctuations dominate. 

In this letter we present a fluctuating LBE (FLBE). This offers a fully consistent dis- 
cretization of the equations of fluctuating nonlinear hydrodynamics for an isothermal fluid, 
opening the way to more accurate and efhcient simulation of many of the mesoscale physics 
problems mentioned previously, such as colloid hydrodynamics. Its generalization to multi- 
component fluids is conceptually straightforward; we pursue this elsewhere [13]. Our work also 
raises broader issues for numerical statistical mechanics: how best to implement fluctuation- 
dissipation theorems (FDTs), derived in the continuum with respect for appropriate conser- 
vation laws, in a system discretized in space and time [14]. We contend that accuracy and 
efficiency are best combined if FDT is made to hold directly on the discretized dynamics. 
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We start with the LBE for a single-component fluid, which can be viewed as a discretization 
of the Boltzmann equation for the coUisional dynamics of a dilute gas [15]: 

/,(x + c„ i + 1) = /,(x, t) + (/,(x, t) - /°(x, t)) (1) 

Here fi represents the local mass density of particles in a phase space cell (x,Ci), and is 
normalized so that X]r=i fi ~ tinid mass density at x. The spatial coordinates x 

are discretized onto a unit lattice; a finite set of n velocities Ci is chosen so that in timestep 
At = 1, the resulting 'streaming' displacements Ax^ = c^Ai are lattice vectors. The local 
momentum density and flux are ga = X^i fi'^ia and IIq^ = J2i fiCiaCip] Greek indices indicate 
Cartesian directions. The equilibrium distribution is conditioned by the local values of p, 
ga and Ha/j. The collision operator dj linearly relaxes the local phase-cell densities towards 
this equilibrium. For well chosen lattices and collision operators, the LBE is known to recover 
the isothermal Navier-Stokes equation in the continuum limit at low Mach number [1]. 

The fi in Eq^are ensemble- averaged local distribution functions. Although these depend 
on temperature (indeed, /° for a fluid at rest is Boltzmann distributed [16]), they describe 
only the mean densities in each phase space cell. Accordingly, they neglect the fluctuating 
stress Saf3 in the equations of fluctuating nonlinear hydrodynamics [19] for an isothermal fluid: 

dtp + ^aga = 0, (2a) 
dtga + V;3n„0 = O, (2b) 

^ap ^ gaVp + pSafJ + Vap-yS^ -y^s + Sap ■ (2c) 

In Eqs[21 Va — ga/ P the local fluid velocity, p is the pressure in a quiescent fluid (given here 
by an ideal gas equation of state p — pc^, with Cs the isothermal sound speed), and rjap-ys is 
a tensor of viscosities. The fluctuating stress Sap is a zero-mean Gaussian random variable 
whose variance, for a fluid at temperature T, is fixed by the fluctuation-dissipation theorem 
(FDT) to be (s„^(x, t)s^s[^ , t')) = 2kBTT^ap^s5{-^ ~ ^')S{t ~ <')• 

An important precursor to our work is that of Ladd [8] . His method consists of adding a 
stochastic piece to the microscopic stress tensor in Eq^ Following through to the continuum 
limit, Eq[21is recovered. However, this process ensures only that FDT holds in the hydrody- 
namic limit k —^ 0. In practice, no numerical algorithm is ever used in this limit, which would 
require infinite computational resources. Indeed, in the noiseless LBE, acceptable hydrody- 
namic behavior of the fluid is maintained right up to /c ~ 2 [20] ; accordingly accurate colloid 
hydrodynamics is achieved with colloids of rather small radius, e.g., R — 2.5 [8]. (Lengths 
are expressed in lattice units, so that fc = 2 has a wavelength of tt). But unless the correct 
noise behavior is implemented over a similarly wide k range, breakdown of FDT at high k 
is liable to infect the whole simulation, once noise is added. We show below, for a simple 
benchmark problem of colloids in traps, that this can indeed occur, resulting in errors of order 
10%. (These errors might be reduced by making R several times larger at fixed volume frac- 
tion; but this is not computationally efficient.) Moreover, as found below, such errors show 
nontrivial dependence on parameters such as the colloid volume fraction. Without exhaustive 
testing, therefore, one cannot know whether Ladd's algorithm is accurate or not, in any given 
region of parameter space. In contrast, our own algorithm appears to give good equilibration 
of colloidal degrees of freedom whenever the suspending fluid is accurately in equilibrium. 

To make progress, we note that in addition to the hydrodynamic degrees of freedom ap- 
pearing in Eq|21 the LBE necessarily involves the dynamics of local, non-hydrodynamic modes 
often called 'ghosts' [21,22]. These are needed in Eq^to maintain isotropic and Gallilean 



^R. Adhikari, ^K. Stratford, ^M.E. Gates and ^A.J. Wagner: Fluctuating lattice BoltzmannS 



invariant hydrodynamics [15,23]. As shown below, the method of Ref. [8] effectively sets T = 
for the ghost mode noise. There is then a risk that ghosts drain thermal fluctuations away 
from the hydrodynamic modes, which therefore may never reach equilibrium. 

To create our improved FLBE, we promote Eq^into a discrete Langevin equation, where 
the fi are interpreted as instantaneous, fluctuating densities in phase space: 

/,(x + c„ t + 1) = /,(x, t) + r y (/,(x, t) - /°(x, i)) + (3) 

with noise terms ^i(x, t) that give fluctuations in the populations in each phase space cell. To 
recover thermal equilibrium, the must be linked, by an FDT, to all sources of dissipation in 
the collision operator Cij. The required FDT must allow that, at a given site and timestep, 
the are correlated in such a way as to exactly conserve p and ga ■ To recover continuum hy- 
drodynamics, the collision process necessarily avoids relaxing the conserved quantities, which 
propagate only via the 'streaming step' x ^ x + c^ in Eqn][l]. Below we diagonalize the coUi- 
sional dissipation, using established methods [23] . In contrast to previous work, we determine 
the accompanying noise structure consistent with statistical mechanical principles. 

First, note that each unnormalised eigenvector m° of £ijwith eigenvalue (— 1/ra) is asso- 
ciated with a corresponding local density M°(x,i) via 

n 

M''(x,t)=^m«/,(x,t),a = l...,i (4) 

i=l 

For example, the density p{'x.,t) equates to A/^(x, t) where m\ = IVi. For a general LBE 
in d dimensions containing n velocities (a 'DdQn model'), there are precisely n eigenvectors, 
corresponding to the n degrees of freedom contained in the fi at a given site. Each 
relaxes at a rate given by the eigenvalue — I/tq. Conserved hydrodynamic variables (mass, 
momentum, but not stress) do not relax; we set the corresponding eigenvalues to zero (—1/ra = 
0) without loss of generality. Note that our method allows for a multi-relaxation time kernel 
but does not assume it; ghost noise is important even with a single relaxation time. 

A complete mode count then consists of the null eigenvector nn\ corresponding to the 
density p; the d null eigenvectors Cia corresponding to the d components of the momentum 
ga] \d{d-{- 1) eigenvectors Qiap = CiaCip — c^Sap corresponding to (independent components 
of) the deviatoric momentum flux Sap = — i^ksTSap, where v is the number density; 
and the remaining n — {1 + d + ^d{d + 1)) ghost mode eigenvectors. The latter are model- 
dependent, but readily computed from Cij in any given implementation of the LBE [1]. Our 
n eigenvectors form a complete, orthogonal basis in the finite velocity space of the LBE. Thus 

/,(x,i) = ^u;,m°M°(x,i)A^" (5) 

a 

with normalizers N°- obeying iV" J^i Wimfnii = Sab and Wi a set of known weights [1]. 

A similar expansion can be applied to any function defined on the velocity space, including 
the noise ^i (x, t) . Let us first set £,i = + with H the hydrodynamic subspace and G is 
its complement, the ghost subspace. Identifying the H eigenvectors as above, we write 

ef = E -.<r(x, t)N^ = Lip ^ + 9i^\ (6) 

where the noise terms associated with p, g^, Sap are denoted by p...; repeated Greek indices 
are summed on. 
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Within this framework, conservation laws have a very simple effect: mass and momentum 
conservation demand p — and = in Eq|3 We are left with 



Here the first term is : it produces thermal fluctuations in the stress tensor, and is the 
noise used by Ladd [8]. The remaining terms are ^p: these maintain thermal equilibrium 
for the ghosts, and are new to this work. Our approach differs crucially from earlier work on 
continuous and discrete stochastic versions of the Boltzmann equation [17,18]; these attempted 
to derive the noise statistics from the collision kernel alone, without taking proper account of 
mass and momentum conservation. 

It remains to determine the covariance matrix of the n — d — 1 nonzero noises (^'^^^) that 
now appear in Eas l6l7l Projecting EqOonto the basis of eigenvectors, taking deviations SM'^ 
around a spatially uniform state and Fourier transforming in space we obtain 



where r°^(k) = exp(-jk • c,) and r'' = {1 - t^^). Eq|Hl represents a set 

of coupled Langevin equations in discrete time. We now assume that all correlations are k- 
independent; this is justified below. Next we square both sides of the above equation, average 
over the noise, and invoke stationarity of equal-time correlators to obtain the required FDT: 



The T-dependence is a standard consequence of the discrete time dynamics in Eq|3[l,8]. The 
d+ 1 null modes corresponding to conserved quantities have no dissipation and hence no noise. 
All remaining modes, including ghosts, have both thermal fluctuations and dissipation; for 
consistent dynamics, they must have noise. Note from Eq|Hlthat setting ti, = 1 for ghosts [8], 
does not decouple the hydrodynamic modes from the ghost noise, except strictly at k = 0. 

To complete our calculation of the noise amplitudes, we now quantify the equilibrium 
thermal fluctuations {5M°-5M^). This requires a thermodynamic model for our fluctuating 
fluid. In keeping with the original thinking behind the LBE [1], and also with its practical 
application (low Mach number), we can choose for this model the thermodynamics of an 
ideal gas. The fluctuation matrix {SM°'6M^) is then k independent, as promised above, 
and computable from knowing that all equilibrium phase-cell occupancies in such a gas obey 
Poisson statistics [19]. The required matrix {SM"-SM^) then follows by a change of basis 
from fi to M°, using Eq^ In implementing our ELBE numerically, we instead transform 
in the other direction via Eq[3 This gives from EqEI a set of correlated phase-cell noises 
^j(x, for use in EqEl As shown above, the noises acting on the n different densities fi at 
site X are not independent: they derive from only n — d — 1 underlying independent noises. 
And, as our eigen-analysis of the collision operator makes clear, they arise from the coUisional 
dissipation of the ^d{d + 1) stress modes and the (n — (1 -I- d + ^d{d + 1)) ghost modes. The 
computational overhead of adding ghost noise is slight (~ 10% on run time). 

The consistency of our ELBE can be assessed by measuring numerically the 'equilibration 
ratio' (ER) for fluctuating hydrodynamic quantities. This is the ratio of a measured variance 
to the one required by the Boltzmann distribution [16] at the temperature T chosen for 
the simulation. We set mass, length, and time units so that p = 1 on an unit lattice. In any 
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Fig. 1 - Left panel; Equilibration ratios for the Fourier modes of the density (o), momentum (A), 
and diagonal (□) and off-diagonal (0) stress components in a D3Q15 model with r = 0.75 shown as 
a function of k with k = fc(l, 1, 1). We set ksT = 1/3000. Open symbols: the present work; closed 
symbols: setting = 0. (Compare Fig. 3(d) of [20] where the correct hydrodynamics, without noise, 
is shown to hold up to fc ~ 2.) Right panel: (|(;2,(k)p) in a D2Q9 model with r = 1, plotted as a 
function of 6 for various values of k, when k = fc(cos 9, sin 9) . The radial coordinate is the equilibriation 
ratio with unity as the solid circle. Lower half plane: the present work (note that symbols overlap). 
Upper half plane: setting = 0. Top to bottom at upper right, k = 0.62, 1.25, 1.88, 2.51. 



quiescent fluid the net thermal momentum g in volume AV has variance (g^) — AV pksT [11]. 
For the ideal gas fc^T = p/v = pc^^/v = (?^/v, where v fixes the phase-cell occupancies referred 
to above, and hence all noise amplitudes. Setting AV = 1, the on-site thermal velocity v = g/ p 
obeys {v^) — ksT. Thus, we must have fc^T ^ = 1/3 (for DdQn models) to satisfy the 
low Mach number requirement of the LBE. 

We now report results for D2Q9 and D3Q15 lattices, using a collision operator with unit 
relaxation time for ghost modes and relaxation time r for stress. The viscosity tensor is that 
of an isotropic fluid, with shear viscosity rj = c^(r — i) and bulk viscosity ( = ^V- (This is 
a particular case of the multirelaxation time operator used in [8].) Fig^ compares ERs for 
different hydrodynamic modes on a D3Q15 lattice, and shows the dependence on direction 
in wavevector space of the momentum fluctuations on a D2Q9 lattice. We also show the 
results found by setting = in Eq[3 as per Ref. [8]. As anticipated above, the latter gives 
acceptable ERs only for fc ^ 1, with values of 0.2 — 0.7 for 1 < k < 2, although the noiseless 
hydrodynamics remains accurate here [20]. In contrast, our algorithm is accurate throughout. 
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Fig. 2 - Equilibration ratios for on site density (□), momentum (o) and stress (A diagonal, <) off- 
diagonal) fluctuations as a function of the relaxation time r in a D3Q15 model. We set fcsT = 1/3000. 
Upper datasets: the present work (note that symbols overlap), lower datasets: setting = 0. 
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Fig. 3 - Equilibration ratios for colloidal displacements in traps as a function of periodic box size 
L for ?7 = 1/6 (left) rj = 1/60 (right) and two colloidal radii (circles R — 1.5, squares R — 2.3) in a 
D3Q15 model. We set ksT = 1/19200. Upper datasets: the present work; lower datasets, the same 
but setting — 0. Error bars are of order the symbol size. 



Figure El shows on-site ERs for various fluctuating quantities as a function of r. These 
involve integration over all k, and could be viewed as a 'worst case' scenario for the local 
temperature seen by a small colloidal particle. Our algorithm is accurate, even in this limit. 
We recommend choosing ksT < 1/3000 so that the fluid is accurately incompressible; we then 
find equilibration at the one percent level or better. 

To simulate colloids directly, we have combined our noise implementation with the al- 
gorithm of [24] for the colloidal particles [25]. We find results for velocity autocorrelators 
that offer a slight improvement over the case with ghost noise omitted [26] ; but for standard 
parameter settings both methods give acceptably accurate values for the colloidal diffusion 
constant D — ksT / {Q-kt^R). A more demanding test of equilibration is to simulate colloids in 
harmonic traps, of the kind often encountered in experiments with optical tweezers [27]. We 
placed N particles, each in a separate confining potential, at regular intervals on an lattice. 
The particles interact hydrodynamically, but this does not affect the Boltzmann distribution 
for their thermal displacements within (well-separated) traps. 

A statistically decisive comparison can be made for N — 1^ when a given particle is 
surrounded not by other colloids in independent traps but by its own periodic images, with 
which it interacts hydrodynamically. In figure |21 we show ERs for the thermal displacement 
of such a particle. Whereas our results are satisfactory, Ladd's algorithm shows worsening 
ERs as L/R is decreased [28]. Any attempt to circumvent these large (> 10%) systematic 
errors arising from omission of ghost noise by defining an 'effective temperature' would merely 
move those errors into D instead. FDT violations in the ghost sector seemingly present major 
obstacles to the accurate simulation of colloids, particularly at modest interparticle spacings 
(modest L/R) [24] which are overcome by our algorithm [29]. 

We have also performed tests with > 1 (now fixing L = 32). The ER is found for each 
particle separately by running for 10^ timesteps with rj — 1/60 and ksT — 1/3000. (The 
equilibration time within a trap is of order r^/D ^ 10'*.) For A^ = 64 particles of i? = 2.3 in 
traps where the correct rms thermal displacement is tq = -\/3, we find ERs of 0.861 ± 0.018 
for Ladd's algorithm and 0.957 ± 0.017 for our own. For R = 2.3 and tq = V3/2 {N = 64) 
the results are 0.928 ± 0.007 versus our 1.055 ± 0.011. For R = 6.23 and ro = \^/2 (A^ = 16) 
the results are 0.962 ±0.011 versus our 1.020 ±0.015 [30]. These findings with iV > 1 suggest 
broadly improved results with our algorithm. However a complete exploration of parameter 
space would be needed to fully resolve the differences of the two methods: without this, one 
cannot know in advance when it is safe to neglect ghost noise. This situation is made worse for 
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dense colloids without traps; here the confining potential created on a colloid by its neighbours 
might not be sampled correctly, but that cannot be checked by explicit calculation. 

In conclusion, by writing a discrete Langevin equation at lattice level and carefully applying 
a fluctuation dissipation theorem, we have derived a fluctuating lattice Boltzmann equation in 
which the ghost modes are fully thermalized. This removes their fluctuation-draining effects 
on the hydrodynamic sector at finite k. Our method gives improved numerical equilibration of 
fluctuating quantities at all k, thus resolving a potentially major obstacle to the use of lattice 
Boltzmann methods in the simulation of thermal fluids, including colloid hydrodynamics. 
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